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We present a numerical study on an interacting monomer-dimer model with nearest neigh- 
bor repulsion on a square lattice, which possesses two symmetric absorbing states. The model 
is observed to exhibit two nearby continuous transitions; the Z2 symmetry-breaking order- 
disorder transition and the absorbing transition with directed percolation criticality. We 
find that the symmetry-breaking transition shows a non-Ising critical behavior, and that the 
absorbing phase becomes critical, in the sense that the critical decay of the dimer density 
observed at the absorbing transition persists even within the absorbing phase. Our findings 
call for further studies on the microscopic models and corresponding continuum description 
belonging to the generalized voter university class. 

I. INTRODUCTION 

Interacting systems often exhibit phase transitions into permanently arrested, absorbing states in 
non-equilibrium (NEQ) steady states [1, 2]. Intense efforts have been devoted to classifying these 
NEQ phase transitions into distinct universality classes, which is one of the most fundamental issues 
in the NEQ statistical mechanics [3]. Although the guiding principles for the classification are still 
lacking in general, symmetry is undoubtedly one of the crucial ingredients in the classification. In 
the most prominent universality class knovirn as the directed percolation (DP) class [4-6], systems 
possess either a single absorbing state or multiple ones not connected by any underlying symmetry 
[7, 8]. Another newly established class known as the generalized voter (GV) universality class [9, 10] 
encompasses a large variety of models possessing two symmetric absorbing states, including statistical 
models for the social dynamics such as opinion spreading [11]. This universality class embraces most 
of one dimensional models belonging to the so called directed Ising class, such as NEQ kinetic Ising 
model [12, 13], interacting monomer-dimer (IMD) model [14], and the voter model with spin-exchange 
[9]. 

A fundamental and unique feature of the GV class is that, due to the symmetry of the absorbing 
states, the models in the class actually involve two types of phase transitions, namely, the Z2 symme- 
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try breaking (SB) order-disorder and the active-absorbing phase transitions. In one dimension, the 
coincidence of these two transitions appears to be inevitable since the SB should be induced by the 
absorbing transition. In two dimensions, the models belonging to the GV class arc shown to exhibit 
two types of behaviors as follows. The SB order-disorder transition and the active-absorbing transition 
can take place simultaneously at a single critical point [15], where the defect density shows a logarith- 
mic decay, as in the original voter model [16]. Close but separate occurrence of these two transitions 
is also observed [17]. Continuum field theoretic descriptions unifying these two types of behavior were 
formulated via appropriate Langevin equations for the local order parameter field [10, 18, 19]. 

We here present an extensive simulation study on the IMD model with nearest neighbor (NN) 
repulsion on a square lattice. We observe that the two transitions are occurring successively at close 
but clearly distinct points. Moreover, we find that the SB transition reveals a new non-Ising critical 
behavior, and that the absorbing transition exhibits a DP-critical behavior (due to the Z2-SB), that 
is, the density of active sites, p{t), shows a power-law decay in time as p{t) ~ t~'^ with <p ~ 0.45, in 
agreement with the corresponding DP exponent in two dimensions. Interestingly, we further find that 
the entire absorbing phase becomes critical in the sense that the same power-law relaxation of p{t) 
persists even deep in the absorbing phase. In previous works for the microscopic models belonging 
to the GV class, we feel that the nature of the SB transition was presumed to be of pure Ising type, 
and thus was not carefully studied so far: the Binder cumulants associated with the SB transition 
took the values U ~ 0.56 [18] or ?7 ~ 0.59 [17], smaller than the universal Ising value Vising — 0.611 
[20]. In contrast, the pure Ising nature of the SB transition in the corresponding continuum Langevin 
description was reported though explicit results are not shown [10]. Our findings thus call for further 
studies on these microscopic models belonging to the GV universality class, and on the corresponding 
Langevin equations. Further related discussions are given below. 

II. INTERACTING MONOMER-DIMER MODEL ON A SQUARE LATTICE 

Here we consider the IMD model with NN repulsion between the same type of particles, which 
is an interesting variant of the original Ziff-Gulari-Barshad's catalytic surface-reaction model [21]. A 
monomer (denoted by A) adsorbs on a (randomly selected) vacant site with free adsorption-attempt 
probability P. Likewise a dimcr (denoted by BB) adsorbs on two (randomly selected) NN vacant 
sites with free adsorption-attempt probability (1 — P)- The parameter P is the single parameter in 
the present model. Adsorbed dimer is assumed to instantly dissociate into two monomers. In two 
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FIG. 1: The Binder cuniulants for the present niodel((a)) and the EQ Ising niodel((b)): U{Psb) = 0.592(1) is 
smaller than the universal value Uismg — 0.611 for the EQ two dimensional Ising model. 

dimensions, one has to carefully implement the kinetic rules respecting the constraint of NN repulsion: 
(a) Monomer A adsorbs on a vacant site only if its NN sites are completely empty or have one or 
more B's regardless of the presence of A. Adsorbed A then reacts with one of the NN B's to form 
the 'molecule' (AB) which leaves the system, (b) Likewise, dimer BB adsorbs on two vacant NN sites 
only if their six NN sites satisfy one of the following conditions: (1) they are completely vacant. (2) 
they have one or more A's in the three left-half (LH) NN sites, and at the same time have no B in 
the three right-half (RH) NN sites, and the vice versa. (3) they have one or more A's on both the LH 
and RH NN sites, regardless of the presence of B. The system thus has the two symmetric absorbing 
states: 'anti- ferromagnetic (AF)' arrangements of vacant sites and A's. 
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FIG. 2: (a) The staggered magnetization Ms{P) vs P for various system sizes L near the SB transition, (b) 
Ms vs L at the SB transition. It shows a power law decrease of Ms{Psb) with L as AIs{Psb) ~ L^^^'^ with 
P/i' = 0.161(3). Ising value fi/v = 0.125 (blue line) is shown for comparison. 



III. SIMULATION RESULTS AND DISCUSSIONS 



Simulation is performed on a square lattice with linear size L using periodic boundary condition. 
Throughout the work, simulation always starts from the vacuum state. The 'spin' value is assigned at 
each site i: = 1, —1,0 for A, vacancy, and B, respectively. Then one can define the local 'staggered 
magnetization' at site i as mj = (— l)*^"'"*«(Tj where ix(y) is the x{y) coordinate of the site i. The order 
parameter for the SB transition is given by the staggered magnetization M^, the average of the local 
staggered magnetization, i.e., Mg = {Ms), where Ms = J2i tni/N {N = L?'). In order to locate the SB 
transition point, we measure the corresponding Binder cumulant U defined asU = \ — {M^) / (3(A^s)'^). 

Shown in Fig. 1(a) is C/ as a function of P for various system sizes L = 64 ~ 256. It takes the 
size-independent universal value U{Psb) = 0.592(1) at the transition, from which the critical point 
is determined to be Psb = 0.6385(1). For comparison, we show U for the EQ two-dimensional Ising 
model in Fig. 1(b) which gives the universal value Uismg — 0.611 at = 2.2692. We see that the 
observed value U{Psb) for the present model is found to be smaller than the pure Ising value. We 
contend that this difference is a manifestation of the non-Ising nature of the SB transition. In order 
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FIG. 3: Finite size scaling results of L'^/^Mg vs L^/^iP - Psb) with ^ = 0.161 and = l-0((a)), and with the 
Ising values /3 — 0.125 and i' = l((b)). 

to investigate this point further, we measure the critical exponents of the SB transition via the finite 
size scaling (FSS) method and a NEQ kinetic method involving the critical coarsening. 

Figure 2(a) shows a plot of Mg vs P for various system sizes. It is expected that Mg obeys the FSS 
of the form M^(P) = L'^/'' F{{P - Psb)L^^'')- Thus its critical value Ms{Psb) will vanish wih system 
size L exhibiting a power law decay as Ms{Psb) ~ L~^^'^. Shown in Fig. 2(b) is a double-log plot of 
Ms{Psb) vs L whose slope for L = 128,160,200,256 gives -13 /u = -0.161(3) (Note that the Ising 
slope —ji/v = —0.125 considerably deviates from the actual slope of the data). With this measured 
value of (3/v = 0.161, we find that the best scaling collapse is obtained for u ~ 1.0, as shown in 
Fig. 3(a). In contrast, when we use the exact Ising exponents (3 = 0.125 and = 1, we see that the 
scaling is systematically violated, as shown in Fig. 3(b). Thus the estimated critical exponents are 
given by /3 ~ 0.161 and v ~ 1.0. 

We also measured the fluctuation of the order parameter x = -^[(-^s) ~ (-^s)^]) which is plotted 
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FIG. 4: (a) Fluctuations of the staggered magnetization x vs P for different system size L. (b) x vs L at the 
SB transition. It shows a power law increase with system size L as x ^ L'^^'^ with "f/v — 1.69(2). Ising value 
j/iy — 1.75 (blue line) is shown for comparison. 

against P for different system sizes in Fig. 4(a). It tends to diverge with increasing L, and is expected 
to obey the FSS x = L'^^'^G{L^^'^{P — Psb))- Thus, x{Psb) will diverge with system size as x{Psb) ~ 
L"'^'^ , which is shown in the inset of Fig. 4(b) with ^jv = 1.69(2) (The Ising slope ^jv = 1.75 is shown 
to deviate from the actual slope of the data). With this measured value, figure 5 (a) shows that the 
best scaling is again obtained for ~ 1.0. Note again that the scaling breaks down for the Ising value 
7 = 1.75 and = 1, as demonstrated in Fig. 5(b). Therefore, we get 7 = 1.69(2) and v ~ 1.0. 

To further check the non-Ising nature of the SB transition, we obtain the critical exponent r/ from 
the dynamic order-parameter correlation function C(r, t) = { J2i 'iT^i{t)^i+r{t)) /N measured during the 
critical coarsening process ensued by a quench from an initial vacuum state to the SB critical point in 
a large system with L = 2048, which does not involve the data collapse of FSS method. Figure 6 shows 
C(r, t) for various times. At long times, the steady state behavior emerges and grows in the short and 
intermediate length regions, exhibiting a critical decay as C(r) ~ r~'^ with rj = 0.310(3). Note that 
the Ising slope rjising = 0.25 shown in the figure considerably deviates from the actual data. We also 
note that the above estimates of the critical exponents are in good agreement with the known scaling 
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FIG. 5: Finite size scaling results of L ^l^x vs L^/^iP ~ Psb) with 7 = 1.69 and ly ~ 1.0((a)), and with the 
Ising values ^jv — 1.75 and v = l((b)). 

relations 2/?/z^ = r/ (in two dimensions) and ^jv = 2 — r] since = 0.322(6) and 7/2^ = 1.69(2). 

Having confirmed tlie non-Ising criticality of tlie SB transition, we now discuss a possible under- 
lying origin for this non-Ising behavior. Although we do not have a solid theoretical argument for 
the observed non-Ising behavior, the AF Ising-like domain morphology appears to provide a crucial 
qualitative hint for the non-Ising nature of the SB transition (see Fig. 7). In the present model, the 
two symmetric absorbing states correspond to the checkerboard pattern of the monomer A's due to 
the NN repulsion and NEQ nature of dynamic rules. One distinct feature in the present model is the 
presence of the dimers. We can see that the existence of the dimers induces (via dynamic rules) longer 
ranged effective repulsion interactions between the monomers, creating wide vacant regions percolat- 
ing through the system (see Fig. 7). From the viewpoint of the AF order, these vacant regions contain 
densely packed defects, and must contain many B elements as well within them. These defects as well 
as B's are expected to divergingly fluctuate in their numbers at the SB transition. Quite probably 
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FIG. 6: The staggered magnetization correlation function C{r,t) vs r for different times t = 10 ■ 2"(n = 
7, 8, • • • , 13) measured during the critical coarsening process {P = Psb)- The steady-state part shows a critical 
decay r~'' with t] = 0.31 (see the inset). Power law decay with Ising value rj = 0.25 (blue line) is shown for 
comparison. Inset: Power law decay of C(r, t) vs r over steady state region for the two latest times, t — 40960 
(diamond) and t = 81920 (circle). The exponent 77 = 0.31 is obtained from the slope measured over the distance 
up to r = 10. 

it is the existence of these defects and B elements (with their diverging fluctuations) that causes the 
non-Ising critical behavior of the SB transition. For example, it is rather clear that these defects and 
B's can play a crucial role such that the spatial correlation of the AF order parameter exhibits more 
rapid algebraic decay compared to the pure Ising case at the SB transition. 

Our preliminary results for the number fluctuation of the defects indeed indicates that it diverges 
at the SB transition, and its presentation is deferred to a later publication. Here we present only the 
number fluctuations of the B element, i.e., xb = ^[{Pb) — (P-b)^] at the SB transition. We indeed 
find that this quantity also exhibits a diverging behavior at the SB transition, though much weaker 
one than the staggered magnetization fluctuations X) a-s shown in Fig. 8. The finite size scaling of xb 
is given by xb = L^^^'^^ H{L^/'^^ {P — Psb))- At the SB transition, the double-log plot xb vs L gives 
the slope 'Jb/^b = 0.212(3). With this value, we find that the best scaling collapse is obtained for 
l/uB = 0.80(5), as shown in Fig. 8. This gives jb = 0.265(22) and ub = 1.25(8). 
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FIG. 7: A typical steady-state configuration near the SB transition (P = 0.6383) in a system with L — 128. 
Blue site denotes monomer state (A), red one dimer state (B), and empty one vacant state. 

In the present model, the SB and the absorbing transitions are split, and the former occurs within 
the active phase (as shown in Fig. 9), and the full AF order is achieved only within the absorbing 
phase. By monitoring the relaxation of the dimer density pBsii), the order parameter of the absorbing 
transition, for various P's, as shown in Fig. 9, we obtain an estimate for the absorbing transition point 
as Pa = 0.6495(10) where the dimer density pBB{t) exhibits a power law decay pBB{t) ~ t~'^ with 
(j) ~ 0.45. This exponent is the same as that in the two-dimensional DP absorbing transition, which 
indicates that the absorbing transition occurring at Pa belongs to the DP universality class. This 
behavior is in accord with the expectation [7, 8] that the the critical behavior is of DP class if the 
Z2 symmetry of the absorbing state is broken. Remarkably, the same critical relaxation of the dimer 
density is observed to persist even within the entire absorbing phase {P > Pa)- In this sense, the 
absorbing phase may be considered as critical. Note that since Z2 symmetry is already broken for 
P > Pa, the system evolves toward the absorbing steady-state via coarsening of the AF domains. The 
AF domains are absorbing and thus monomers and dimers cannot fall on them. They therefore must 
fall on the vacant regions which A's and B's are scattered within. Coarsening process thus is driven by 
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FIG. 8: Lower inset: The dimer fluctuations xb vs P for different system size L. Upper inset: xb vs L at tfie 
SB transition. It shows a power law increase with system size L as xb U^bI^b ^b/vb — 0.212(3). Main 
frame: FSS of L^^^/^^xs vs L'^/''"{P - Psb) with 7b = 0.265 and vb = 1-25. 



this interface-alone-fluctuations, and is expected to exhibit a dynamic scaling for the spatio-temporal 
correlations of the AF order parameter. We may consider the observed persistent critical decay of the 
dimer density even within the absorbing phase as associated with this inevitable coarsening dynamics 
from initial vacuum state evolving toward the final (symmetry-broken) absorbing state at all points of 
the absorbing phase. In other words, if we assume that the total number of the dimers is proportional 
to the volume of the vacant regions, then the power law decay of the dimers will ensue from the power 
law decay of the volume of vacant regions associated with the kinetic self-similarity of the system 
during the coarsening (i.e., phase ordering) process. As for the universal DP exponent value in the 
whole absorbing phase, we may draw some analogy to the "zero temperature fixed point" in the phase 
ordering of a standard Ising model quenched to a temperature below the transition temperature Tc, 
where the phenomenology of coarsening at any temperature below Tc can be essentially represented 
by the same universality of dynamic criticality. 

We now discuss on the continuum stochastic description for the GV universality class. In Ref.[10], 
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FIG. 9: The dimer density relaxation pBB{t) vs t for various values of P in a system with L = 2048. The 
absorbing transition point is estimated to be Pa = 0.6495(10), at which pBsit) exhibits a power law relaxation 
PBsit) ^ t^'^ with (p ~ 0.45 (dotted line denotes a line with the slope -0.45). The same power-law relaxation is 
observed to persist within the absorbing phase. 

the following Langevin equation for the local SB order parameter field <j) was proposed for a unified 
description of the two types of behaviors alluded earlier for the models in the GV class: 

dt<p = (ac^ - b(j)^){l - (t?) + V^<t) + a^l^i (III.l) 

where ^ is a uncorrelated Gaussian noise with zero mean and unit variance, and a, 6, and a are 
constants. This equation is a modified version of the Langevin equation for the voter model [22, 23] 
((III.l) with a = 6 = 0) so as to incorporate the splitting of the Z2 SB transition and the absorbing 
transition, which was first observed in the model studied in [17]. Since the order parameter field 
fluctuates around = (and hence (1 — cfP') ~ 1) near the SB transition, the multiplicative factor in 
the noise term in (III.l) disappears, reducing (III.l) to that of the model A. This argument seems to 
be consistent with the report of Ref. [10] that the convincing numerical solution of (III.l) reveals that 
the Z2 SB transition described by (III.l) belongs to the Ising universality class, which are evidenced 
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by the FSS scaling and measurements of the Binder cumulant and time auto correlation function for 
the order parameter (though these results are not shown exphcitly in that paper). It is thus quite 
probable that the Langevin equation (III.l) cannot fully describe the non-Ising critical behavior of the 
SB transition observed in the present model. 

Recently, Vazquez and Lopez (VL) [18] derived a similar Langevin equation from the microscopic 
model of the Ising spins Gr on a square lattice with spin flip dynamics controlled by a flip probability 
which is given by a function /(^) of the local magnetization ip{v,t) = ^r)/Z {Z being the 

coordination number) with Yl^, denoting sum over Z neighboring spins of Uj-. This was done by 
mapping the master equation into an approximate Fokker-Planck equation, and in turn into the 
corresponding Langevin equation which is given by 

X(0) = [a + c+{d-2a-Zb)4)\ 

Y{(^) = {l-4>'^){c + d4>'^) + {a-c + 2d)4>V^4> (IIL2) 

The corresponding spin-flip probability function is given by /('0) = {l + il)){c+ail) + dilP' — hil)^)/2 with 
a, 6, c and d being constants. VL found through the simulation on a square lattice that the splitting 
of the SB and absorbing transitions is observed when the next-next-nearest-neighbor [Z = 12) spins 
are taken into account in the calculation of the local magnetization ip, and that the Binder cumulant 
U ~ 0.56 was obtained at the SB transition. We tend to believe that this value of U, considerably 
smaller than the pure Ising value U = 0.61, indicates that the SB transition exhibits a non-Ising critical 
behavior. It is interesting to see that unlike (III.l), the corresponding Langevin equation (III. 2) does 
not seem to reduce to the model A near the SB transition due to the forms of X{(f)) and Y{<p) in (III.2). 
This feature could leave a room for the possibility that the Langevin equation (III.2) might be able 
to generate a non-Ising critical behavior of the SB transition, which deserves a thorough investigation 
on the critical behavior of the SB transitions in both microscopic model and corresponding Langevin 
equation proposed by VL. 

IV. SUMMARY 

In the present work, we investigated a two dimensional IMD model with short-range repulsion, 
which has two symmetric absorbing states. We find that the SB and absorbing transitions take place 
at two nearby points, as in some models [17, 18] with extended interaction range belonging to the GV 
universality class. Interestingly, we observe numerically that the SB transition exhibits a new non-Ising 
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critical behavior with the exponents (3 ~ 0.161, 7 ~ 1.69, 2± 1.0, and rj ~ 0.31. We also find that 
the absorbing transition is of the DP-type, with a remarkable feature that the dimer density shows 
the same critical decay (with the DP-exponent) within the entire absorbing phase. The non-Ising 
critical behavior of the SB transition is considered to be closely related to the fact that the effective 
extended range of repulsion between monomers mediated by dimers (via NEQ dynamic rules) creates 
wide vacant regions between the AF domains, which can be considered as closely packed defects from 
the view of the AF order. We suggested that the diverging fluctuations of these defects and dimers can 
lead to a non-Ising critical behavior of the AF order parameter at the SB transition. The persistent 
criticality in the absorbing phase can be attributed to the self-similar nature of the coarsening process 
(i.e., the AF domain growth) evolving toward symmetry-broken absorbing state. Finally, we raise a 
possibility that a recently derived Langcvin equation based on the spin-flip dynamics of a collection of 
Ising spins on a square lattice with a general flip probability might be able to describe similar non-Ising 
critical behavior of the SB transition. It would thus be valuable to carry out further analysis on those 
models and continuum stochastic equations showing separate transitions in the GV universality class. 
It would be also interesting to conduct more detailed investigations into NEQ kinetics in both active 
and absorbing phases in the present model. 
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